package beast.evolution.substitutionmodel;

import beast.core.Description;
import beast.evolution.datatype.Aminoacid;
import beast.evolution.datatype.DataType;
import beast.math.statistic.DiscreteStatistics;

@Description("WAG model of amino acid evolution by " +
        "S. Whelan and N. Goldman. 2001. Molecular biology and evolution 18.5 (2001): 691-699 ")

public class WAG extends EmpiricalSubstitutionModel {

    @Override
    double[][] getEmpiricalRates() {
        double[][] rate = new double[20][20];

        // Q matrix from Beast 1
        rate[0][1] = 0.610810;
        rate[0][2] = 0.569079;
        rate[0][3] = 0.821500;
        rate[0][4] = 1.141050;
        rate[0][5] = 1.011980;
        rate[0][6] = 1.756410;
        rate[0][7] = 1.572160;
        rate[0][8] = 0.354813;
        rate[0][9] = 0.219023;
        rate[0][10] = 0.443935;
        rate[0][11] = 1.005440;
        rate[0][12] = 0.989475;
        rate[0][13] = 0.233492;
        rate[0][14] = 1.594890;
        rate[0][15] = 3.733380;
        rate[0][16] = 2.349220;
        rate[0][17] = 0.125227;
        rate[0][18] = 0.268987;
        rate[0][19] = 2.221870;

        rate[1][2] = 0.711690;
        rate[1][3] = 0.165074;
        rate[1][4] = 0.585809;
        rate[1][5] = 3.360330;
        rate[1][6] = 0.488649;
        rate[1][7] = 0.650469;
        rate[1][8] = 2.362040;
        rate[1][9] = 0.206722;
        rate[1][10] = 0.551450;
        rate[1][11] = 5.925170;
        rate[1][12] = 0.758446;
        rate[1][13] = 0.116821;
        rate[1][14] = 0.753467;
        rate[1][15] = 1.357640;
        rate[1][16] = 0.613776;
        rate[1][17] = 1.294610;
        rate[1][18] = 0.423612;
        rate[1][19] = 0.280336;

        rate[2][3] = 6.013660;
        rate[2][4] = 0.296524;
        rate[2][5] = 1.716740;
        rate[2][6] = 1.056790;
        rate[2][7] = 1.253910;
        rate[2][8] = 4.378930;
        rate[2][9] = 0.615636;
        rate[2][10] = 0.147156;
        rate[2][11] = 3.334390;
        rate[2][12] = 0.224747;
        rate[2][13] = 0.110793;
        rate[2][14] = 0.217538;
        rate[2][15] = 4.394450;
        rate[2][16] = 2.257930;
        rate[2][17] = 0.078463;
        rate[2][18] = 1.208560;
        rate[2][19] = 0.221176;

        rate[3][4] = 0.033379;
        rate[3][5] = 0.691268;
        rate[3][6] = 6.833400;
        rate[3][7] = 0.961142;
        rate[3][8] = 1.032910;
        rate[3][9] = 0.043523;
        rate[3][10] = 0.093930;
        rate[3][11] = 0.533362;
        rate[3][12] = 0.116813;
        rate[3][13] = 0.052004;
        rate[3][14] = 0.472601;
        rate[3][15] = 1.192810;
        rate[3][16] = 0.417372;
        rate[3][17] = 0.146348;
        rate[3][18] = 0.363243;
        rate[3][19] = 0.169417;

        rate[4][5] = 0.109261;
        rate[4][6] = 0.023920;
        rate[4][7] = 0.341086;
        rate[4][8] = 0.275403;
        rate[4][9] = 0.189890;
        rate[4][10] = 0.428414;
        rate[4][11] = 0.083649;
        rate[4][12] = 0.437393;
        rate[4][13] = 0.441300;
        rate[4][14] = 0.122303;
        rate[4][15] = 1.560590;
        rate[4][16] = 0.570186;
        rate[4][17] = 0.795736;
        rate[4][18] = 0.604634;
        rate[4][19] = 1.114570;

        rate[5][6] = 6.048790;
        rate[5][7] = 0.366510;
        rate[5][8] = 4.749460;
        rate[5][9] = 0.131046;
        rate[5][10] = 0.964886;
        rate[5][11] = 4.308310;
        rate[5][12] = 1.705070;
        rate[5][13] = 0.110744;
        rate[5][14] = 1.036370;
        rate[5][15] = 1.141210;
        rate[5][16] = 0.954144;
        rate[5][17] = 0.243615;
        rate[5][18] = 0.252457;
        rate[5][19] = 0.333890;

        rate[6][7] = 0.630832;
        rate[6][8] = 0.635025;
        rate[6][9] = 0.141320;
        rate[6][10] = 0.172579;
        rate[6][11] = 2.867580;
        rate[6][12] = 0.353912;
        rate[6][13] = 0.092310;
        rate[6][14] = 0.755791;
        rate[6][15] = 0.782467;
        rate[6][16] = 0.914814;
        rate[6][17] = 0.172682;
        rate[6][18] = 0.217549;
        rate[6][19] = 0.655045;

        rate[7][8] = 0.276379;
        rate[7][9] = 0.034151;
        rate[7][10] = 0.068651;
        rate[7][11] = 0.415992;
        rate[7][12] = 0.194220;
        rate[7][13] = 0.055288;
        rate[7][14] = 0.273149;
        rate[7][15] = 1.486700;
        rate[7][16] = 0.251477;
        rate[7][17] = 0.374321;
        rate[7][18] = 0.114187;
        rate[7][19] = 0.209108;

        rate[8][9] = 0.152215;
        rate[8][10] = 0.555096;
        rate[8][11] = 0.992083;
        rate[8][12] = 0.450867;
        rate[8][13] = 0.756080;
        rate[8][14] = 0.771387;
        rate[8][15] = 0.822459;
        rate[8][16] = 0.525511;
        rate[8][17] = 0.289998;
        rate[8][18] = 4.290350;
        rate[8][19] = 0.131869;

        rate[9][10] = 3.517820;
        rate[9][11] = 0.360574;
        rate[9][12] = 4.714220;
        rate[9][13] = 1.177640;
        rate[9][14] = 0.111502;
        rate[9][15] = 0.353443;
        rate[9][16] = 1.615050;
        rate[9][17] = 0.234326;
        rate[9][18] = 0.468951;
        rate[9][19] = 8.659740;

        rate[10][11] = 0.287583;
        rate[10][12] = 5.375250;
        rate[10][13] = 2.348200;
        rate[10][14] = 0.462018;
        rate[10][15] = 0.382421;
        rate[10][16] = 0.364222;
        rate[10][17] = 0.740259;
        rate[10][18] = 0.443205;
        rate[10][19] = 1.997370;

        rate[11][12] = 1.032220;
        rate[11][13] = 0.098843;
        rate[11][14] = 0.619503;
        rate[11][15] = 1.073780;
        rate[11][16] = 1.537920;
        rate[11][17] = 0.152232;
        rate[11][18] = 0.147411;
        rate[11][19] = 0.342012;

        rate[12][13] = 1.320870;
        rate[12][14] = 0.194864;
        rate[12][15] = 0.556353;
        rate[12][16] = 1.681970;
        rate[12][17] = 0.570369;
        rate[12][18] = 0.473810;
        rate[12][19] = 2.282020;

        rate[13][14] = 0.179896;
        rate[13][15] = 0.606814;
        rate[13][16] = 0.191467;
        rate[13][17] = 1.699780;
        rate[13][18] = 7.154480;
        rate[13][19] = 0.725096;

        rate[14][15] = 1.786490;
        rate[14][16] = 0.885349;
        rate[14][17] = 0.156619;
        rate[14][18] = 0.239607;
        rate[14][19] = 0.351250;

        rate[15][16] = 4.847130;
        rate[15][17] = 0.578784;
        rate[15][18] = 0.872519;
        rate[15][19] = 0.258861;

        rate[16][17] = 0.126678;
        rate[16][18] = 0.325490;
        rate[16][19] = 1.547670;

        rate[17][18] = 2.763540;
        rate[17][19] = 0.409817;

        rate[18][19] = 0.347826;

        // Current (May 2011) version from http://www.ebi.ac.uk/goldman/WAG/index.html
        //
//			0.551571 
//			0.509848  0.635346 
//			0.738998  0.147304  5.429420 
//			1.027040  0.528191  0.265256  0.0302949 
//			0.908598  3.035500  1.543640  0.616783  0.0988179 
//			1.582850  0.439157  0.947198  6.174160  0.021352  5.469470 
//			1.416720  0.584665  1.125560  0.865584  0.306674  0.330052  0.567717 
//			0.316954  2.137150  3.956290  0.930676  0.248972  4.294110  0.570025  0.249410 
//			0.193335  0.186979  0.554236  0.039437  0.170135  0.113917  0.127395  0.0304501 0.138190 
//			0.397915  0.497671  0.131528  0.0848047 0.384287  0.869489  0.154263  0.0613037 0.499462  3.170970 
//			0.906265  5.351420  3.012010  0.479855  0.0740339 3.894900  2.584430  0.373558  0.890432  0.323832  0.257555 
//			0.893496  0.683162  0.198221  0.103754  0.390482  1.545260  0.315124  0.174100  0.404141  4.257460  4.854020  0.934276 
//			0.210494  0.102711  0.0961621 0.0467304 0.398020  0.0999208 0.0811339 0.049931  0.679371  1.059470  2.115170  0.088836  1.190630 
//			1.438550  0.679489  0.195081  0.423984  0.109404  0.933372  0.682355  0.243570  0.696198  0.0999288 0.415844  0.556896  0.171329  0.161444 
//			3.370790  1.224190  3.974230  1.071760  1.407660  1.028870  0.704939  1.341820  0.740169  0.319440  0.344739  0.967130  0.493905  0.545931  1.613280 
//			2.121110  0.554413  2.030060  0.374866  0.512984  0.857928  0.822765  0.225833  0.473307  1.458160  0.326622  1.386980  1.516120  0.171903  0.795384  4.378020 
//			0.113133  1.163920  0.0719167 0.129767  0.717070  0.215737  0.156557  0.336983  0.262569  0.212483  0.665309  0.137505  0.515706  1.529640  0.139405  0.523742  0.110864 
//			0.240735  0.381533  1.086000  0.325711  0.543833  0.227710  0.196303  0.103604  3.873440  0.420170  0.398618  0.133264  0.428437  6.454280  0.216046  0.786993  0.291148  2.485390 
//			2.006010  0.251849  0.196246  0.152335  1.002140  0.301281  0.588731  0.187247  0.118358  7.821300  1.800340  0.305434  2.058450  0.649892  0.314887  0.232739  1.388230  0.365369  0.314730 
        //
//			0.0866279 0.043972  0.0390894 0.0570451 0.0193078 0.0367281 0.0580589 0.0832518 0.0244313 0.048466  0.086209  0.0620286 0.0195027 0.0384319 0.0457631 0.0695179 0.0610127 0.0143859 0.0352742 0.0708956
        //
        //
//			 A   R   N   D   C   Q   E   G   H   I   L   K   M   F   P   S   T   W   Y   V
//			Ala Arg Asn Asp Cys Gln Glu Gly His Ile Leu Lys Met Phe Pro Ser Thr Trp Tyr Val
        //
//			  
//				Symmetrical part of the WAG rate matrix and aa frequencies,
//			estimated from 3905 globular protein amino acid sequences forming 182
//			protein families.
//				The first part above indicates the symmetric 'exchangeability'
//			parameters, where s_ij = s_ji.  The s_ij above are not scaled, but the
//			PAML package will perform this scaling.
//				The second part gives the amino acid frequencies (pi_i)
//			estimated from the 3905 sequences.  The net replacement rate from i to
//			j is Q_ij = s_ij*pi_j.
//				Prepared by Simon Whelan and Nick Goldman, December 2000.
        //
//			Citation:
//			Whelan, S. and N. Goldman.  2001.  A general empirical model of
//				protein evolution derived from multiple protein families using
//				a maximum likelihood approach.  Molecular Biology and
//				Evolution 18:691-699.

/*
		// Q matrix
		rate[0][1] = 0.551571; rate[0][2] = 0.509848; 
		rate[0][3] = 0.738998; rate[0][4] = 1.027040; 
		rate[0][5] = 0.908598; rate[0][6] = 1.582850; 
		rate[0][7] = 1.416720; rate[0][8] = 0.316954; 
		rate[0][9] = 0.193335; rate[0][10] = 0.397915; 
		rate[0][11] = 0.906265; rate[0][12] = 0.893496; 
		rate[0][13] = 0.210494; rate[0][14] = 1.438550; 
		rate[0][15] = 3.370790; rate[0][16] = 2.121110; 
		rate[0][17] = 0.113133; rate[0][18] = 0.240735; 
		rate[0][19] = 2.006010; 

		rate[0][1] = 0.551571;rate[0][2] = 0.509848;
		rate[0][3] = 0.738998;rate[0][4] = 1.027040;
		rate[0][5] = 0.908598;rate[0][6] = 1.582850;
		rate[0][7] = 1.416720;rate[0][8] = 0.316954;
		rate[0][9] = 0.193335;rate[0][10] = 0.397915;
		rate[0][11] = 0.906265;rate[0][12] = 0.893496;
		rate[0][13] = 0.210494;rate[0][14] = 1.438550;
		rate[0][15] = 3.370790;rate[0][16] = 2.121110;
		rate[0][17] = 0.113133;rate[0][18] = 0.240735;
		rate[0][19] = 2.006010;

		rate[1][2] = 0.635346;rate[1][3] = 0.147304;
		rate[1][4] = 0.528191;rate[1][5] = 3.0355;
		rate[1][6] = 0.439157;rate[1][7] = 0.584665;
		rate[1][8] = 2.13715;rate[1][9] = 0.186979;
		rate[1][10] = 0.497671;rate[1][11] = 5.35142;
		rate[1][12] = 0.683162;rate[1][13] = 0.102711;
		rate[1][14] = 0.679489;rate[1][15] = 1.22419;
		rate[1][16] = 0.554413;rate[1][17] = 1.16392;
		rate[1][18] = 0.381533;rate[1][19] = 0.251849;


		rate[2][3] = 5.42942;rate[2][4] = 0.265256;
		rate[2][5] = 1.54364;rate[2][6] = 0.947198;
		rate[2][7] = 1.12556;rate[2][8] = 3.95629;
		rate[2][9] = 0.554236;rate[2][10] = 0.131528;
		rate[2][11] = 3.01201;rate[2][12] = 0.198221;
		rate[2][13] = 0.0961621;rate[2][14] = 0.195081;
		rate[2][15] = 3.97423;rate[2][16] = 2.03006;
		rate[2][17] = 0.0719167;rate[2][18] = 1.086;
		rate[2][19] = 0.196246;

		rate[3][4] = 0.0302949;rate[3][5] = 0.616783;
		rate[3][6] = 6.17416;rate[3][7] = 0.865584;
		rate[3][8] = 0.930676;rate[3][9] = 0.039437;
		rate[3][10] = 0.0848047;rate[3][11] = 0.479855;
		rate[3][12] = 0.103754;rate[3][13] = 0.0467304;
		rate[3][14] = 0.423984;rate[3][15] = 1.07176;
		rate[3][16] = 0.374866;rate[3][17] = 0.129767;
		rate[3][18] = 0.325711;rate[3][19] = 0.152335;


		rate[4][5] = 0.0988179;rate[4][6] = 0.021352;
		rate[4][7] = 0.306674;rate[4][8] = 0.248972;
		rate[4][9] = 0.170135;rate[4][10] = 0.384287;
		rate[4][11] = 0.0740339;rate[4][12] = 0.390482;
		rate[4][13] = 0.39802;rate[4][14] = 0.109404;
		rate[4][15] = 1.40766;rate[4][16] = 0.512984;
		rate[4][17] = 0.71707;rate[4][18] = 0.543833;
		rate[4][19] = 1.00214;

		rate[5][6] = 5.46947;rate[5][7] = 0.330052;
		rate[5][8] = 4.29411;rate[5][9] = 0.113917;
		rate[5][10] = 0.869489;rate[5][11] = 3.8949;
		rate[5][12] = 1.54526;rate[5][13] = 0.0999208;
		rate[5][14] = 0.933372;rate[5][15] = 1.02887;
		rate[5][16] = 0.857928;rate[5][17] = 0.215737;
		rate[5][18] = 0.22771;rate[5][19] = 0.301281;


		rate[6][7] = 0.567717;rate[6][8] = 0.570025;
		rate[6][9] = 0.127395;rate[6][10] = 0.154263;
		rate[6][11] = 2.58443;rate[6][12] = 0.315124;
		rate[6][13] = 0.0811339;rate[6][14] = 0.682355;
		rate[6][15] = 0.704939;rate[6][16] = 0.822765;
		rate[6][17] = 0.156557;rate[6][18] = 0.196303;
		rate[6][19] = 0.588731;

		rate[7][8] = 0.24941;rate[7][9] = 0.0304501;
		rate[7][10] = 0.0613037;rate[7][11] = 0.373558;
		rate[7][12] = 0.1741;rate[7][13] = 0.049931;
		rate[7][14] = 0.24357;rate[7][15] = 1.34182;
		rate[7][16] = 0.225833;rate[7][17] = 0.336983;
		rate[7][18] = 0.103604;rate[7][19] = 0.187247;


		rate[8][9] = 0.13819;rate[8][10] = 0.499462;
		rate[8][11] = 0.890432;rate[8][12] = 0.404141;
		rate[8][13] = 0.679371;rate[8][14] = 0.696198;
		rate[8][15] = 0.740169;rate[8][16] = 0.473307;
		rate[8][17] = 0.262569;rate[8][18] = 3.87344;
		rate[8][19] = 0.118358;

		rate[9][10] = 3.17097;rate[9][11] = 0.323832;
		rate[9][12] = 4.25746;rate[9][13] = 1.05947;
		rate[9][14] = 0.0999288;rate[9][15] = 0.31944;
		rate[9][16] = 1.45816;rate[9][17] = 0.212483;
		rate[9][18] = 0.42017;rate[9][19] = 7.8213;


		rate[10][11] = 0.257555;rate[10][12] = 4.85402;
		rate[10][13] = 2.11517;rate[10][14] = 0.415844;
		rate[10][15] = 0.344739;rate[10][16] = 0.326622;
		rate[10][17] = 0.665309;rate[10][18] = 0.398618;
		rate[10][19] = 1.80034;

		rate[11][12] = 0.934276;rate[11][13] = 0.088836;
		rate[11][14] = 0.556896;rate[11][15] = 0.96713;
		rate[11][16] = 1.38698;rate[11][17] = 0.137505;
		rate[11][18] = 0.133264;rate[11][19] = 0.305434;


		rate[12][13] = 1.19063;rate[12][14] = 0.171329;
		rate[12][15] = 0.493905;rate[12][16] = 1.51612;
		rate[12][17] = 0.515706;rate[12][18] = 0.428437;
		rate[12][19] = 2.05845;

		rate[13][14] = 0.161444;rate[13][15] = 0.545931;
		rate[13][16] = 0.171903;rate[13][17] = 1.52964;
		rate[13][18] = 6.45428;rate[13][19] = 0.649892;


		rate[14][15] = 1.61328;rate[14][16] = 0.795384;
		rate[14][17] = 0.139405;rate[14][18] = 0.216046;
		rate[14][19] = 0.314887;

		rate[15][16] = 4.37802;rate[15][17] = 0.523742;
		rate[15][18] = 0.786993;rate[15][19] = 0.232739;


		rate[16][17] = 0.110864;rate[16][18] = 0.291148;
		rate[16][19] = 1.38823;

		rate[17][18] = 2.48539;rate[17][19] = 0.365369;


		rate[18][19] = 0.31473;
*/


////        WAG* matrix, a variant of the WAG matrix from the same paper
////		
////		0.589718 
////		0.514347  0.67416 
////		0.731152  0.159054  5.30821 
////		1.21324   0.568449  0.233527   0.0379056 
////		1.03344   3.02808   1.62299    0.657364   0.0999068 
////		1.55788   0.443685  1.00122    6.04299    0.0284956  5.6037 
////		1.41993   0.629768  1.12717    0.88357    0.312544   0.346823  0.588609 
////		0.317684  2.31211   3.9337     0.958529   0.341479   4.87366   0.599188   0.279542 
////		0.214596  0.187262  0.527321   0.0390513  0.198958   0.125999  0.124553   0.0310522  0.162975 
////		0.400822  0.51821   0.144354   0.0869637  0.451124   0.873266  0.154936   0.067443   0.508952  3.1554 
////		0.881639  5.74119   2.88102    0.480308   0.0719929  4.19125   2.45392    0.381514   0.854485  0.320597  0.255092 
////		0.887458  0.660816  0.198404   0.0992829  0.428648   1.64018   0.294481   0.184545   0.40117   3.94646   4.81956   0.877057 
////		0.213179  0.122792  0.0848492  0.0458258  0.485001   0.109241  0.0873936  0.0552962  0.631713  1.06458   2.10414   0.0832422  1.14516 
////		1.51861   0.711498  0.204905   0.444152   0.109081   0.913179  0.720567   0.254626   0.722123  0.111722  0.422851  0.588203   0.179858  0.165205 
////		3.52499   1.35611   3.90127    1.09965    1.35221    0.87908   0.822025   1.33618    0.876688  0.321774  0.351913  1.05314    0.554077  0.563999  1.54694 
////		2.24161   0.594177  2.06787    0.395176   0.522957   0.829315  0.889765   0.236489   0.54992   1.48876   0.351564  1.45173    1.56873   0.188237  0.802531  4.02507 
////		0.135395  1.24086   0.0746093  0.142159   0.728065   0.208163  0.176397   0.366467   0.261223  0.259584  0.706082  0.159261   0.565299  1.58681   0.135024  0.528249  0.118584 
////		0.270321  0.386714  1.05269    0.326191   0.481954   0.210494  0.209621   0.108982   4.31772   0.44009   0.427718  0.155623   0.437069  6.49269   0.212945  0.742154  0.286443  2.42261 
////		1.92496   0.282892  0.193323   0.155419   1.10899    0.32893   0.588443   0.190095   0.119749  7.48376   1.82105   0.300343   2.03324   0.653015  0.325745  0.23769   1.4088    0.396884  0.353358 
////
////		0.0866279 0.043972  0.0390894 0.0570451 0.0193078 0.0367281 0.0580589 0.0832518 0.0244313 0.048466  0.086209  0.0620286 0.0195027 0.0384319 0.0457631 0.0695179 0.0610127 0.0143859 0.0352742 0.0708956
////
////
////		 A   R   N   D   C   Q   E   G   H   I   L   K   M   F   P   S   T   W   Y   V
////		Ala Arg Asn Asp Cys Gln Glu Gly His Ile Leu Lys Met Phe Pro Ser Thr Trp Tyr Val
////
////		  
////			Symmetrical part of the WAG* rate matrix and aa frequencies,
////		estimated from 3905 globular protein amino acid sequences forming 182
////		protein families.
////			The first part above indicates the symmetric 'exchangeability'
////		parameters, where s_ij = s_ji.  The s_ij above are not scaled, but the
////		PAML package will perform this scaling.
////			The second part gives the amino acid frequencies (pi_i)
////		estimated from the 3905 sequences.  The net replacement rate from i to
////		j is Q_ij = s_ij*pi_j.
////			Prepared by Simon Whelan and Nick Goldman, December 2000.
////
////		Citation:
////		Whelan, S. and N. Goldman.  2001.  A general empirical model of
////			protein evolution derived from multiple protein families using
////			a maximum likelihood approach.  Molecular Biology and
////			Evolution 18:691-699.		
//
//		rate[0][1] = 0.589718;rate[0][2] = 0.514347;
//		rate[0][3] = 0.731152;rate[0][4] = 1.21324;
//		rate[0][5] = 1.03344;rate[0][6] = 1.55788;
//		rate[0][7] = 1.41993;rate[0][8] = 0.317684;
//		rate[0][9] = 0.214596;rate[0][10] = 0.400822;
//		rate[0][11] = 0.881639;rate[0][12] = 0.887458;
//		rate[0][13] = 0.213179;rate[0][14] = 1.51861;
//		rate[0][15] = 3.52499;rate[0][16] = 2.24161;
//		rate[0][17] = 0.135395;rate[0][18] = 0.270321;
//		rate[0][19] = 1.92496;
//
//		rate[1][2] = 0.67416;rate[1][3] = 0.159054;
//		rate[1][4] = 0.568449;rate[1][5] = 3.02808;
//		rate[1][6] = 0.443685;rate[1][7] = 0.629768;
//		rate[1][8] = 2.31211;rate[1][9] = 0.187262;
//		rate[1][10] = 0.51821;rate[1][11] = 5.74119;
//		rate[1][12] = 0.660816;rate[1][13] = 0.122792;
//		rate[1][14] = 0.711498;rate[1][15] = 1.35611;
//		rate[1][16] = 0.594177;rate[1][17] = 1.24086;
//		rate[1][18] = 0.386714;rate[1][19] = 0.282892;
//
//
//		rate[2][3] = 5.30821;rate[2][4] = 0.233527;
//		rate[2][5] = 1.62299;rate[2][6] = 1.00122;
//		rate[2][7] = 1.12717;rate[2][8] = 3.9337;
//		rate[2][9] = 0.527321;rate[2][10] = 0.144354;
//		rate[2][11] = 2.88102;rate[2][12] = 0.198404;
//		rate[2][13] = 0.0848492;rate[2][14] = 0.204905;
//		rate[2][15] = 3.90127;rate[2][16] = 2.06787;
//		rate[2][17] = 0.0746093;rate[2][18] = 1.05269;
//		rate[2][19] = 0.193323;
//
//		rate[3][4] = 0.0379056;rate[3][5] = 0.657364;
//		rate[3][6] = 6.04299;rate[3][7] = 0.88357;
//		rate[3][8] = 0.958529;rate[3][9] = 0.0390513;
//		rate[3][10] = 0.0869637;rate[3][11] = 0.480308;
//		rate[3][12] = 0.0992829;rate[3][13] = 0.0458258;
//		rate[3][14] = 0.444152;rate[3][15] = 1.09965;
//		rate[3][16] = 0.395176;rate[3][17] = 0.142159;
//		rate[3][18] = 0.326191;rate[3][19] = 0.155419;
//
//
//		rate[4][5] = 0.0999068;rate[4][6] = 0.0284956;
//		rate[4][7] = 0.312544;rate[4][8] = 0.341479;
//		rate[4][9] = 0.198958;rate[4][10] = 0.451124;
//		rate[4][11] = 0.0719929;rate[4][12] = 0.428648;
//		rate[4][13] = 0.485001;rate[4][14] = 0.109081;
//		rate[4][15] = 1.35221;rate[4][16] = 0.522957;
//		rate[4][17] = 0.728065;rate[4][18] = 0.481954;
//		rate[4][19] = 1.10899;
//
//		rate[5][6] = 5.6037;rate[5][7] = 0.346823;
//		rate[5][8] = 4.87366;rate[5][9] = 0.125999;
//		rate[5][10] = 0.873266;rate[5][11] = 4.19125;
//		rate[5][12] = 1.64018;rate[5][13] = 0.109241;
//		rate[5][14] = 0.913179;rate[5][15] = 0.87908;
//		rate[5][16] = 0.829315;rate[5][17] = 0.208163;
//		rate[5][18] = 0.210494;rate[5][19] = 0.32893;
//
//
//		rate[6][7] = 0.588609;rate[6][8] = 0.599188;
//		rate[6][9] = 0.124553;rate[6][10] = 0.154936;
//		rate[6][11] = 2.45392;rate[6][12] = 0.294481;
//		rate[6][13] = 0.0873936;rate[6][14] = 0.720567;
//		rate[6][15] = 0.822025;rate[6][16] = 0.889765;
//		rate[6][17] = 0.176397;rate[6][18] = 0.209621;
//		rate[6][19] = 0.588443;
//
//		rate[7][8] = 0.279542;rate[7][9] = 0.0310522;
//		rate[7][10] = 0.067443;rate[7][11] = 0.381514;
//		rate[7][12] = 0.184545;rate[7][13] = 0.0552962;
//		rate[7][14] = 0.254626;rate[7][15] = 1.33618;
//		rate[7][16] = 0.236489;rate[7][17] = 0.366467;
//		rate[7][18] = 0.108982;rate[7][19] = 0.190095;
//
//
//		rate[8][9] = 0.162975;rate[8][10] = 0.508952;
//		rate[8][11] = 0.854485;rate[8][12] = 0.40117;
//		rate[8][13] = 0.631713;rate[8][14] = 0.722123;
//		rate[8][15] = 0.876688;rate[8][16] = 0.54992;
//		rate[8][17] = 0.261223;rate[8][18] = 4.31772;
//		rate[8][19] = 0.119749;
//
//		rate[9][10] = 3.1554;rate[9][11] = 0.320597;
//		rate[9][12] = 3.94646;rate[9][13] = 1.06458;
//		rate[9][14] = 0.111722;rate[9][15] = 0.321774;
//		rate[9][16] = 1.48876;rate[9][17] = 0.259584;
//		rate[9][18] = 0.44009;rate[9][19] = 7.48376;
//
//
//		rate[10][11] = 0.255092;rate[10][12] = 4.81956;
//		rate[10][13] = 2.10414;rate[10][14] = 0.422851;
//		rate[10][15] = 0.351913;rate[10][16] = 0.351564;
//		rate[10][17] = 0.706082;rate[10][18] = 0.427718;
//		rate[10][19] = 1.82105;
//
//		rate[11][12] = 0.877057;rate[11][13] = 0.0832422;
//		rate[11][14] = 0.588203;rate[11][15] = 1.05314;
//		rate[11][16] = 1.45173;rate[11][17] = 0.159261;
//		rate[11][18] = 0.155623;rate[11][19] = 0.300343;
//
//
//		rate[12][13] = 1.14516;rate[12][14] = 0.179858;
//		rate[12][15] = 0.554077;rate[12][16] = 1.56873;
//		rate[12][17] = 0.565299;rate[12][18] = 0.437069;
//		rate[12][19] = 2.03324;
//
//		rate[13][14] = 0.165205;rate[13][15] = 0.563999;
//		rate[13][16] = 0.188237;rate[13][17] = 1.58681;
//		rate[13][18] = 6.49269;rate[13][19] = 0.653015;
//
//
//		rate[14][15] = 1.54694;rate[14][16] = 0.802531;
//		rate[14][17] = 0.135024;rate[14][18] = 0.212945;
//		rate[14][19] = 0.325745;
//
//		rate[15][16] = 4.02507;rate[15][17] = 0.528249;
//		rate[15][18] = 0.742154;rate[15][19] = 0.23769;
//
//
//		rate[16][17] = 0.118584;rate[16][18] = 0.286443;
//		rate[16][19] = 1.4088;
//
//		rate[17][18] = 2.42261;rate[17][19] = 0.396884;
//
//
//		rate[18][19] = 0.353358;
//		

        return rate;
    }

    @Override
    public double[] getEmpiricalFrequencies() {
        double[] f = new double[20];
        // frequencies from May 2011
        f[0] = 0.0866279;
        f[1] = 0.043972;
        f[2] = 0.0390894;
        f[3] = 0.0570451;
        f[4] = 0.0193078;
        f[5] = 0.0367281;
        f[6] = 0.0580589;
        f[7] = 0.0832518;
        f[8] = 0.0244313;
        f[9] = 0.048466;
        f[10] = 0.086209;
        f[11] = 0.0620286;
        f[12] = 0.0195027;
        f[13] = 0.0384319;
        f[14] = 0.0457631;
        f[15] = 0.0695179;
        f[16] = 0.0610127;
        f[17] = 0.0143859;
        f[18] = 0.0352742;
        f[19] = 0.0708956;
        // frequencies from Beast 1
        f[0] = 0.0866;
        f[1] = 0.0440;
        f[2] = 0.0391;
        f[3] = 0.0570;
        f[4] = 0.0193;
        f[5] = 0.0367;
        f[6] = 0.0581;
        f[7] = 0.0833;
        f[8] = 0.0244;
        f[9] = 0.0485;
        f[10] = 0.0862;
        f[11] = 0.0620;
        f[12] = 0.0195;
        f[13] = 0.0384;
        f[14] = 0.0458;
        f[15] = 0.0695;
        f[16] = 0.0610;
        f[17] = 0.0144;
        f[18] = 0.0353;
        f[19] = 0.0709;
        return f;
    }

    @Override
    public int[] getEncodingOrder() {
        Aminoacid dataType = new Aminoacid();
        String codeMap = dataType.getCodeMap();
        int[] codeMapNrs = new int[dataType.getStateCount()];
        String encoding = "ARNDCQEGHILKMFPSTWYV";
        for (int i = 0; i < dataType.getStateCount(); i++) {
            codeMapNrs[i] = encoding.indexOf(codeMap.charAt(i));
        }
        return codeMapNrs;
    }

    @Override
    public boolean canHandleDataType(DataType dataType) {
        return dataType instanceof Aminoacid;
    }

    public static void main(String[] args) {

        WAG wag = new WAG();

        //String aminoAcids = "ARNDCQEGHILKMFPSTWYV";
        boolean[] class1 = {false,true,false,false,true,true,true,false,false,true,true,
                false,true,false,false,false,false,true,true,true};

        int within1 = 0;
        int within2 = 0;
        int between = 0;
        double[] w1 = new double[45];
        double[] w2 = new double[45];
        double[] b = new double[100];

        double[][] rates = wag.getEmpiricalRates();

        for (int i = 0; i < 20; i++) {
            for (int j = i+1; j < 20; j++) {
                if (class1[i]) {
                    if (class1[j]) {
                        w1[within1] = rates[i][j];
                        within1 += 1;
                    } else {
                        b[between] = rates[i][j];
                        between += 1;
                    }
                } else {
                    if (!class1[j]) {
                        w2[within2] = rates[i][j];
                        within2 += 1;
                    } else {
                        b[between] = rates[i][j];
                        between += 1;
                    }
                }
            }
        }

        System.out.println("Within 1 mean rate = " + DiscreteStatistics.mean(w1));
        System.out.println("Within 2 mean rate = " + DiscreteStatistics.mean(w2));
        System.out.println("Between mean rate = " + DiscreteStatistics.mean(b));

        System.out.println("Within 1 rate stdev = " + DiscreteStatistics.stdev(w1));
        System.out.println("Within 2 rate stdev = " + DiscreteStatistics.stdev(w2));
        System.out.println("Between rate stdev = " + DiscreteStatistics.stdev(b));

        System.out.println("Within 1 rate stderr = " + DiscreteStatistics.stdev(w1)/Math.sqrt(within1));
        System.out.println("Within 2 rate stderr = " + DiscreteStatistics.stdev(w2)/Math.sqrt(within2));
        System.out.println("Between rate stderr = " + DiscreteStatistics.stdev(b)/Math.sqrt(between));

        double sse = 0;
        double meanb = DiscreteStatistics.mean(b);
        for (int i = 0; i < b.length; i++) {
            sse += (meanb - b[i]) * (meanb - b[i]);
        }
        sse /= b.length;
        sse = Math.sqrt(sse);

        System.out.println("Between rate sse = " + sse);

        System.out.println("Within 1 count = " + within1);
        System.out.println("Within 2 count = " + within2);
        System.out.println("Between count = " +  between);
    }



} // class WAG
